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ABSTRACT 



The rTpplication of Mitrovic's Method enables the control systems 
engineer to analyze and synthesize control systems as a function of two 
parameters. The feedback compensation of sampled-data control systems 
is thoroughly investigated as a multi-parameter problem, with emphasis 
placed on the design of the compensation in the Mitrovic parameter plane. 
The minimization of the control system settling time is shown to be parti- 
cularly adapted to the parameter plane. Several examples are included 
which demonstrates this usefulness of the parameter plane in solving the 
minimum settling time problem. A stability relationship is devised for 
digital computer determination of dynamic system stability which permits 
the formulation of an n-dimensioned parameter space. Applications of this 
digital creation of an n-dimensioned space are considered. This ability 
to design feedback compensation for the sampled-data system as a function 
of more than one parameter is shoxTn to lead to greater accuracy in closed 
loop root positioning; and in certain cases, to permit the design con- 
sideration of several constraints simultaneously. 
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CHAPT^TR T 



GEmm]. CONCKPTS AND THEORY 

^ ^ Int: roduc t: ion , 

The design of multivariable feedback compensation for a feedback 
control system, linear or sampled-data, has heretofore been accomplished 
largely by one of two methods. Either repeated calculations using the 
classical one parameter design techniques, such as root locus, has been 
used; or, if a more accurate study is required, the designer could resort 
to analog simulation and a *'trial and error'’ approach to the design pro- 
blem. Dusan Mitrovic, {3], in his early work with the algebraic methods 
for locating the roots of a polynomial introduced a new analysis technique 
that permitted analysis of the root location in terras of two variable co- 
efficients. Later, Siljak, [4] and [5], devised additional algebraic 
methods whereby two parameters could be selected from a portion of any, 
or all, of the polynomial coefficients, and the root analysis performed 
in terms of these parameters. This extension has proven exceptionally 
useful in tv;o parameter analysis of control systems whose characteristic 
equations have coefficients that are each functions of the same two vari- 
ables or system parameters. 

The object of this thesis, then is to apply the results of Mitrovic 
and Siljak to the design of feedback compensation for the sampled-data 
control system. A secondary objective is to show the superiority of the 
two parameter method over the more classical techniques, and to demon- 
strate the adaptability of the entire design problem to solution by digital 
computer. 
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feed- 



Let us consider the characteristic equation of the sampled-data 
back control system as. 



or, 



F(?.) = 1 + GH(?.) = 0 
^ k 

F(z) = 2 __ ^ 

K*0 



(la) 

(lb) 



where the a *s are constant or parameter varying coefficients, and z is 

K 

one of the complex z-plane locations of the closed loop roots of the system. 
Now> consider the generalized point in the z-plane> 

z . (2) 

where lzj= 

and 



j 0 

arg z = e*^ 



= cos ^ - i sin 0 

• ^ - r/ > 






If this general coordinate is substituted into equation (lb), the character- 
istic equation may then be written as: 



F(z) 




+ J O ““ 




(3) 



Equation (3) is, in general, a non-zero quantity and assumes a zero value 
only when the value of and ^ assumed corresponds to that of one of 
the closed loop roots of the system. Thus, if we require that F(z) = 0, 
equation (3) describes the closed loop pole locations if the are 

constants, and every point on the z-plane root locus if one of the 
is parameter varying. Therefore, 



K-. a 

maps the characteristic equation from the z-plane to the Mitrovic plane if 
the a *s are constants, and the parameter plane if the a. 's are variable. 
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Only the plane will be discussed hereafter. 



1 .? ConformaJ_^ Mapping Techniques. 

The requirement that equation (4) be identically equal to zero 
implies that both the real and imaginary parts are zero. Siljak, [5], 
demonstrated that this magnitude requirement on equation (4) leads to 
the two simultaneous equations, 

^ \ \( - 0 <5a) 

K ^ c 

and , 

P «k V S, > - 0 • <5b) 

J\ - t . 

where T^( ^ ) and U^( ) are the Chebyshev functions of the first and 

second kind respectively. Furthermore, since 



V?;) • S’. V ' - \-i< S'.) 

eqxiations (5a) and (5b) can be rewritten as 

l! - \ 1^2 “k-l< S. ) = ” 



( 6 ) 



and 



K-_o 



p \ co>,(c;.) -0 



(7a) 



(7b) 



respectively, with the U, ( vT.)'s generated by the recursive Chebyshev 
formula 



U 



k + ] 






with 



u. ( 1 



u ,(C, ) -1 

- 1 ^Te. 

Siljak, [4] and [5], considered the case in which the are linear 

functions of two scalar parameters, 
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\ ■'^"k + + ‘‘k 



( 9 ) 



The substitution of equation (9) into equations (7a) and (7b) forms the 
set of simultaneous equations. 



and 



where 












(lOa) 



IT . ®2= 



K = o 



K^O 



K*0 



4 - S' ^ <=kto“vs^> 



K«0 



^ -s “k.i<SJ > S' 



K = C? K«0 

Equations (10a) and (10b) can then be solved simultaneously for the two 
parameters, alpha and beta as; 



oC- 






h ^2 ■ ® 2^1 




«2^ 



«l'=2 






«2S 



(11) 



Equations (11) are the fundamental results of Slljak, [41 and [5], 
and provide the basis for the work which follows. Equation (11) can be 
interpreted as the calculated values of the two parameters alpha and beta 
required to put £ closed loop root at a particularly chosen location on 
the z -plane. Thus, by fixing and varying ^ from -1 to +1, upper 

semi-circle contours of constant 6^^ can be mapped on the parameter plane. 
Thus, the argument of each upper semi-circle complex root conformally maps 
onto a single point on the parameter plane, the coordinates of which are 
the values of the two parameters producing that particular root configura- 
tion. Similiarly, in equation (11), can be fixed while C^^is varied 
from zero to infinity producing loci of constant ^ , or angle, on the 



parameter plane. 
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A third type of conformal mapping is yet required to completely specify 
all the roots of the characteristic equation on the parameter plane. The 
loci of constant and constant determine all the complex roots. 

The loci of points on the parameter plane corresponding to constant real 
roots on the z-plane can be obtained directly from the characteristic equa- 
tion, equation (lb). The substitution of equation (9) into equation (lb), 
and letting z = , yields 




which for a given 0^ is a straight line on the parameter plane, regard- 
less of the order of the characteristic equation. 
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CHAPTER II 



SECOND ORDER SYSTEMS 



The second order sampled-data control system provides an excellent 
model for developing the feedback compensation techniques from applica- 
tions of the equations of Chapter I. For the initial design application, 
consider the pure second order system shown in figure (la) . The selected 
design goal is to minimize the settling time by suitable choice of the feed- 
back parameters, a^ and If the states of the system of figure (la) are 

defined as: 



then the system can be represented by the signal flow-graph of figure (lb), 
and the vector-matrix equation 



and 



x^(t) = c(t) 
X2<t) = S(t) 



(13) 



X[(k + 1)T] = 5 X(kT) + A R(kT), 



(14) 



where 



1 



(a^KT^)/2 



T - 



5 



- (a^KT) 



1 



(a2KT) 



(15) 



and 




A 



(16) 



KT 



The characteristic equation of the system is; 




2.1 Constant Magnitude Loci . 



The duration of the transient response, or settling time for a 



linear system of the second order is, [6] 
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Figure 1 

Pure Second Order System 
Die grams 




I 

(a) Block Diagram 



( 17 ) 



For a stable system then, the contours of lines of constant settling time 
of a second order system in the s-plane are 









or s = - 4/t 

s 


• 


(18) 


Accordingly, from one definition of the 2 -transformation, [2], 




Ts 

2 = e 


> 


(19) 



we can determine corresponding contours of constant settling time in the 
2 -plane as. 



-4T/t 

2 = e s 



( 20 ) 



which, since T and t are constants for any particular design, are circles 

s 

-4T/t 

centered at the origin of the 2 -plane of radius e s. 

Thus, the loci of constant 63^ contours on the parameter plane can 

be interpreted as loci of constant settling times; and, when 63^ and T 

are specified, the appropriate t^ can be uniquely determined. Solving 

equation ( 20 ) for t , 

s 



tg = - 4T/ln 6 J 2 . (21) 

Equation (21), however, must be interpreted in terms of the z-plane, dis- 
crete, characteristics. For example, for a step input, the minimum set- 
tling time obtainable is nT seconds, for an n-th order system, when all 
the closed loop roots are at the origin of the z-plane. Accordingly, 
when (Og approaches zero in equation (21), t^ should approach nT. Thus, 
for proper interpretation, equation ( 21 ) must be written as 

t = nT - 4T/ln cO, . (22) 

s ^ 

The particular conclusion of interest here is, in terms of the parameter 
plane, to minimi 2 e t^ one has only to minimi 2 e within the stable 
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region of the plane. 

To demonstrate the effectiveness of the constant 6.^, loci in aid- 
ing the designer meet settling time criteria, consider again the pure 
second order system of figure (1), whose characteristic equation is given 
on page 6. It is first necessary to define the system parameters, so let 

o(. S a^KT^ (23) 

^ S a^KT 

By reference to equation (9) and the characteristic equation, then 




d2 = 2 = 0 C2 = 0 

which when substituted into equation (10) yields 

Bi - -0_j( S,) ‘ S > 

<=1 = <=2 = S .) 

SJ -2 “l( Sj 

“2 - ^\u,( S, ) + 2 

These factors can now be substituted into equation (11) and the constant 
6^2 loci plotted with the aid of a desk calculator or a digital computer. 
Note that if a digital computer is used, as is most likely, a limiting 
process is required to remove the common factor of (0^: ^he numera- 

tor and denominator of the parametric equations when attempting to locate 
the = 0 point on the parameter plane. Figure (2) shows the resulting 

lines of constant on the parameter plane for this system. 

From figure (2), minimum settling time is predicted when 0, 

or when 

= 1.0 , 

and =1.5 
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In order to determine the system response as predicted by equation (14), 
let 

T = K = I.O , 

then = 1.0 

and = 1.5 is required. 

Table 1 lists the calculated magnitude values of the states of the system 
at consecutive sampling intervals using these values of the two parameters. 
The time response of the states is plotted on figure (3). A unit step in- 
put is assumed. 

TABLE 1 



nT X^(t) X 2 <t) 



0 0 0 

1 0.5 1.0 

2 1.0 0 

3 1.0 0 

4 1.0 0 



Since in this example, n = 2, these results are in agreement with that 
predicted by equation (22), and the response is the minimum prototype, 
ripple free response. 

In order to demonstrate the most important design advantages in using 
the parameter plane for feedback compensation; assume that, for some other 
reason, one of the feedback coefficients is constant, say a^ is fixed at 
2.25. Let K = T = 1.0, also be constant. With these constraints on the 
design problem, the parameter plane affords an excellent potential for 
setting the remaining degree of freedom to minimize the settling time, 
in so far as possible v;ith feedback compensation alone. Referring to 
figure (2) for = 2.25, is minimum at 6)^ = 0.5 when yS = 1.875. 

From equation (22), 
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Figure 3 

Response of Pure Second Order System 
with Minimum Settling Time Compensation 
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r= 2T - 4T/1.n 0.5 

or t = 7.8T 

s 

is the predicted settling time. Table 2 lists the state variable magni- 
tudes at the indicated sampling interval in response to a unit step input 
utilizing the above parameter values in equation (14). Figure (4) shows 
the time response of the state variables of the system as a function of 
the sampling interval, or time. 



TABLE 2 



nT 


xpt) 


X2<t) 


0 


0 


0 


\ 


0.5000 


1 . 0000 


2 


0.5000 


-1.0000 


3 


0.3750 


0.7500 


h 


0.5000 


-0.5000 


5 


0.4063 


0.3125 


6 


0.4688 


-0.1875 


7 


0.4297 


0.1094 


8 


0.4531 


-0.0625 


9 


0.4395 


0.0352 


10 


0.4473 


-0.0195 


11 


0.4429 


0.0107 


12 


0.4453 


-0.0059 


13 


0.4440 


0.0032 


14 


0.4447 


-0.0017 


15 


0.4443 


0.0009 


16 


0,4445 


-0.0005 


17 


0.4444 


0.0003 


18 


0.4445 


-0.0001 


19 


0 . 4444 


0.0001 


20 


0.4444 


-0.0000 



The results indicated in Table 2 appear to discredit the validity of 
equation (22) in predicting settling time. Actually, equation (22) simply 
nrovides an engineering approximation to the settling time problem. Analy- 
sis of the results of Table 2 indicates the output position is within 1.957o 
of the final value after 8T seconds, although an actual zero velocity and 
final position of 0.4444 was not achieved until the 41st sampling interval. 
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Figure 

Response of Pure Second Order System 
with Constraints 
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A similar tost was performed with a = 2.25 and a = 1.50 for = 0.8. 
Enue^ition (22) predicts the settling time as t^ = 20T, at which time the 
output had settled to within 3.37o of the final value. In that example, 
zero actual velocity did not occur until t = 99T. However, the settling 
time prediction of equation (22) remains within generally accepted good 
engineering design proximity of the desired value. It must be borne in 
mind ho\>^ever, that for different from 1.0, the conventional definition 
of system error cannot be applied. Fortunately, in many cases, input/out- 
put scaling can be applied which permits the position voltage feedback to 
be set at other than unity. 

2.2 Constant Angle Loci . 

Equation (11), the parametric equations describing 






^ 1*^2 ■ ^ 2^1 






® 2 ^ 




® 2 ^ 






® 1°2 






( 11 ) 



the loci of the complex roots on the parameter plane are functions of both 
and . The constant magnitude curves were obtained by fixing 

and allowing to vary as the mapping parameter. Similiarly, 

can be fixed while varies from zero to infinity to map lines of 
constant z-plane angle for the complex roots onto the parameter plane. 

While the loci of this constant angle on the parameter plane does not have 
any direct usefulness comparable to the loci of constant magnitude and the 
determination of settling time, they do lead to a very important secondary 
property. With loci of constant 63^ and constant drawn on the para- 

meter plane, every point on the z-plane now has a unique point defined on 
the parameter plane in the case of the second order systems. The fact that 
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for higher order systems, several z-plane points may map onto one para- 
meter plane point will be discussed later. 

For the second order case, this one-to-one correspondence between 
points on the z-plane and points on the parameter plane is useful in con- 
formal mapping of various loci from the z-plane to the parameter plane. 

Among the possibilities are the logarithmic constant damping spirals corres- 
ponding to lines of constant damping in the s-plane. Another use would be 
the mapping of lines of constant overshoot, [2] and [7], from the z-plane 
to the parameter plane. This z-plane mapping permits the designer to 
select feedback compensation which simultaneously satisfies overshoot and 
settling time criteria for the second order system. 

2,3 Constant Real Root Loci . 

In order to completely specify all possible z-plane root locations on 
the parameter plane, real root configurations must be allowed for as well 
as the complex root configurations. As indicated by equation (12), the 
loci of constant real roots are straight lines on the parameter plane. 

For the second order system, the closed loop poles are either both real or 
both complex. Therefpre, the zone of real root description and the zone of 
complex root description cannot overlap, except at a critical point where 
the real root is about to emerge from the real axis. 

This zone separation characteristic greatly simplifies the sketching 
of the parameter plane for the second order system. Fortunately, only one 
set of calculations is necessary. The constant 60^ curves must be cal- 
culated. However, since these loci of constant magnitude are straight 
lines, only a minimum amount of calculations is involved. In plotting 
the constant ^ curves, one has only to connect the coordinates of 
equal on the 60,, curves. Then, to plot the lines of constant , 
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( 24 ) 



since from equation (2) we can write, 



0^- S. 



one has only to draw straight lines tangent to the I I == 1*0 line at 
the cO intersections. For example, a line drawn tangent to the ^ = -1.0 
line at the intersection of the =0.5 line describes the real root 



locations 6^ = -0.5. This sketching technique will not necessarily be 
extendable to the higher order systems unfortunately; however, it does 
increase the simplicity of the second order parameter plane. 

Figure (5) is the parameter plane for the pure second order system 
of figure 2 which now has some of the loci of constant angle (solid lines) 
and some of the loci of constant real roots (broken lines) superimposed 
upon it. Note that this parameter plane configuration also accurately 
describes the system stability in terms of the two parameters. 

2.4 Pole Effect . 

In order to study the pole effect of the second order systems on the 
parameter plane, the system shown in figure (6) was studied. This system 
can also be represented by the vector-matrix equation (14), where in this 
case 



p2 ' ' 




i 
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Figure 6 

General Second Order 
System Diagrams 





(b) Signal Plow-Graph 



and 



^,(PT - I + ) 



A 



I- 

The characteristic equation is 





( 25 ) 



where 



= Ka^(T/p-l/p+e’^^/p^) + Ka 2 (l/p-e'P^/p) -1-e 



and 



Aq = Kaj^(l/p^-Te’'’^/p-e ^^/p^) + Ka^Ce *^^/p-l/p)+e 



One important difference between the pure second order system and the 
system with damping becomes evident by inspecting the characteristic equa- 
tion above. In the pure second order case, by judicious definition of the 
system parameters, equation (23), it was possible to produce a parameter 
plane as a function (indirectly) of all physical parameters, K, T, a^, and 
By moving one of the poles into the left-half of the s-plane, it is 
no longer possible to obtain parameters involving T. However, the system 
gain, K, remains a common factor in the characteristic equation coeffici- 
ents. Therefore, for the second order system with damping, we can define 
the parameters as 



For further study of the pole effects, let p = T » 1.0, then from equa- 
tion (9), 




(26) 



d^ « 0.368 
dj^ = -I. 368 
d2 - 1.0 



bp = 0.264 

bj = 0.368 



Cq » -0.632 
Cj = 0.632 



b 



2 



0 
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and from equations (10) , 

= -0.264U_j( - 0.368 ) 

Cj = 0.632U_j( ) Cj » 0.632 ^ ) 

Dj = -0.368U.j( Uj( 

”2- -1-368cO^,Uj(5;) + idlv^(. 5^). 

Figure (7) is the parameter plane for this system, with loci plotted 
for comparison with figure (5). Note that the shifting of one pole from 
the origin of the s~plane into the left-half plane has produced a counter- 
clockwise rotation and a compression of the upper semi-circle contours. 

This suggests that the further into the left half plane of the real pole 
(closer to the origin of the z-plane) the more sensitive the system re- 
sponse is to changes in Additionally, from figure (7), the extremities 

of the complex root region on the parameter plane have been extended, to 
the point where negative values of beta are now within the stable region. 
From figure (7), minimum prototype response to a step input should be ex- 
pected for 

c< = 1.58 
/S = 1.24 

Consequently, there are two ways to achieve minimum prototype compensa- 
tion. First, if K is fixed to some value other than 1.58, for example 
1.0, then non-unity position feedback must be used. With K = 1.0, the 
magnitudes of the state variables at consecutive sampling intervals in 
response to a unit step input, calculated using equation (14), are shown 
by table 3. Note that this response is identical to that shoim on figure 
(3), except for a constant scale factor, which is naturally equal to 
1/a^. 
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TABLE 3 



nT 



X^(t) 



X2(t) 



0 

1 

2 

3 

4 



0 

0.368 

0.632 

0.632 

0.632 



0 

0.632 

0 

0 

0 



For the particular case where the conventional error definition is 
required, to obtain the minimum prototype response, K must be restricted 
to the value of 1.58, Hence, the corresponding value for a 2 is 1.24/K or 
0,785. The calculated response using these parameter settings is shown 
by table 4, assuming a unit step input again. 

TABLE 4 

nT X^(t) X 2 (t) 

0 0 0 

1 0.582 1.0 

2 1.000 0 

3 1.000 0 

4 1.000 0 



2,5 Zero Effect . 

To study the effect of a real zero on the second order parameter 
plane, the system shown in figure (8) was analyzed. For this system, 
the state transition matrices are; 

1 - Kaj^(T + ZT^/2) T - Ka^CT + ZT^/2) 

1 - Ka2<l + ZT) 

K(T + ZT^/2) 

K(1 + ZT) 






and 



- Ka^(l + ZT) 



A 
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Figure 8 

Second Order System 
with Zero, 
Diagrams 





(b) Signal Plow-Graph 



Again, the characteristic equation is 

2 

z + A^z + Aq = 0 , 

where 

Aq = Ka^(ZT^ - 2T) - Ka2(2ZT + 2) + 1 
and Aj = Ka^(ZT^ + 2T) + Ka2(2ZT + 2) - 2 

Again, a choice of parameters for the parameter plane can be made as was 
performed in section 2,4. For further study, Z and T must be fixed in 
magnitude. For this example then, let Z ~ T = lcO» Following the now 
familiar calculations, the parameter plane of figure (9) results. Similiar 
to the results shown for the pole effect, the parameter plane loci have 
been rotated counter-clockwise from the pure second order configuration of 
figure (5). Apparently this counter-clockwise rotation can be associated 
with the relative stability of the system, as both the system with the zero 
and the system with the left-half s-plane pole are inherently more stable 
than the pure second order servo. However, no measure for quantitatively 
describing this relative stability has been discovered. 

From figure (9), minimum prototype, ripple free response to a step in- 
put is predicted for 

C7< = 0.5 

= 0.125 

Again there are two choices for the physical parameters in achieving 
this compensation. If K is variable, unity position feedback may be used. 
If K is fixed, then non-unity position feedback is required. Responses 
similiar to that achieved in the last section result for each of these 
options . 
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CHAPTER III 



HIGHER ORDER SYSTEMS 
3.1 Parameter Limitations . 

It can be readily demonstrated that in order to achieve complete free- 
dom in locating the closed loop roots of an n-th order characteristic 
equation, the control system error signal must contain feedback intelli- 
gence on each of the n dynamically independent **states*^ of the systemo 
Consequently, the system analog or mathematical model consists of n para- 
meters and the parameter space becomes an n-space. Unfortunately, the 
existing Mitrovic, or algebraic, methods have their basis in the solution 
of two simultaneous equations, equations (5a) and (5b), This restricts 
the mathematical treatment of any n-th order system to only two degrees of 
freedom. While this reduction in the degrees of freedom restricts the 
ability to position the closed loop poles in the z-plane, a large class 
of control systems can be reduced to the two parameter problem. 

Typically, the gain setting, K, in figure (10) is dictated by the 
requirements for steady state error. A second degree of freedom is normal- 
ly eliminated by the adoption of the conventional error definition which 
demands unity feedback. Finally, in most physical systems, noise free 
derivative signals higher than the second derivative are unobtainable. 

Thus, the only degrees of freedom for feedback compensation typically avail- 
able are the magnitudes of the first and second derivative feedback. 

With this restriction on the positioning of the closed loop roots in 
mind then, the feedback compensation of systems with order greater than two 
can be designed on the parameter plane using techniques already discussed. 
Two additional features of the parameter plane must now be considered. 
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Typical Block Diagram for 
n-th Order Sampled-Data System 
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Unfortunately, at the present time, these two features tend to reduce 
the effectiveness of the Mitrovic Method. However, it is believed that 
this reduced effectiveness is still an improvement over the existing one 
parameter design techniques. 

3.2 Maximum Modulus Uncertainity . 

The parametric equations describing the parameter plane loci were de- 
rived from the system characteristic equation. Any specified location (M- 
point) on the resulting parameter plane can be interpreted as the values 
of the t\70 parameters which will put a closed loop root at the specified 
position of the z-plane. Therefore, for a typical fourth order system 
for example, the closed loop root configurations shown by figures (11a) and 
(11b) each correspond to a unique point on the parameter plane, each lying 
on the (0., = 1.0 loci. In order to distinguish between the different con- 
figurations on the parameter plane, it is also necessary to plot the 0.5 
and 6^ = 1.5 loci. The intersection of two of the loci of constant 
then adequately describes the location of each pair of complex roots. If 
figure (12) represents the parameter plane for the fourth order system of 
figure (11), then the root location depicted by figure (11a) corresponds 
to point A, while the configuration of figure (11b) corresponds to point 
B on the parameter plane. 

Thus, for an n-th order system, the modulus of each real root or each 
pair of complex roots is completely determined on the parameter plane. 

The accuracy of that determination, however, is directly related to the 
number of loci of constant drawn. In the second order case, it was 

fairly easy to extrapolate between adjacent loci. As the order of the 
characteristic equation increases, the difficulty of making this extra- 
polation increases. It should be noted that this problem of interpreting 
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intersections of constant loci also exists in the second order case when 
the plant gain is insufficient to produce complex roots. Uncertainity 
arises as the order of the system increases. For example, given a fifth 
order system, and a corresponding parameter plane M-point defined by the 
intersection of two constant loci and one real root loci. Obviously, 

this M-point must determine two pair of complex roots and one real root. 
However, with existing techniques, it is impossible to associate any 
specific with a specific complex root. 

3.3 Discontinuities on the Parameter Plane . 

A second feature which may exist in higher order systems involves dis- 
continuities on the constant and loci. Equation (11) describes 

the parametric equations for alpha and beta as a ratio of two polynomials. 
Consequently, whenever 

a discontinuity exists on the loci of constant 63^ or ^ . At the present 

time, the physical significance of this discontinuity is not understood, 
in fact one may not even exist. Its presence, however, does complicate 
the sketching and interpretation of the parameter plane. 

3.4 Dominance and Settling Time . 

In the case of the linear, continuous system, the response of a 
higher order system can often be discussed in terms of a dominant pair of 
roots in the s-plane. This is a pair of complex conjugate roots whose 
residue and/or transient settling time leads to by far the greatest 
portion of the observed physical response. In the s-plane, these roots 
have an argument and a real part considerably less in magnitude than any 
other pair of complex, or real, roots. If such a pair of roots exist, the 
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response can be quite accurately predicted in terms of a second order 
approximation and the settling time predicted by equation (17)- 

Similiarly, for the sampled-data system, the dominance of the system 
response can be discussed in terms of the pair of complex roots, or real 
root, whose modulus or argument is considerably greater than that of any 
other root, real or complex- If such a root exists, the settling time 
can be closely approximated by a slight modification to equation (22) such 
that 

t = nT - 4T/ln ( ) - (27) 

When the dominance cannot be clearly determined by the difference in the 
various arguments, equation (27) becomes less and less accurate in pre- 
dicting the system settline time* 

Alternatively, the minimization of the settling time can be achieved 
qualitatively by minimizing the modulus of each root or pair of complex 
roots. That is to say, compensation must be chosen which will place all 
the z-plane roots within a circle of minimum radius. If this circle of 
minimum radius happens to be the origin of the z-plane, then minimum pro- 
totype response is guaranteed. If, however, all the z-plane roots cannot 
be placed at the origin due to the limitations on the degrees of freedom 
available, minimum settling time will be somewhat greater than minimum 
prototype, and will not be ripple free. On the parameter plane, this 
minimization of the settling time involves the examination of the inter- 
sections of the constant loci for each permissible M-point, and 

selecting that one which does place all the roots within the circle of 
minimum possible radius. Observations of the resulting arguments will 
enable the designer to determine if equation (27) is then valid or not. 
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3 .5 Third Order Example. 



Consider the system shown on figure (13). This system can be re- 
presented by the signal flow graph of figure (14), and the vector-matrix 
equation 



X [(k + 1 )t] = § J((kT) + 4 e(kT) 



where 






1 (3/2 - 2e"'^ + l/2e'^^) (1/2 - + l/2e"^^) 



-T -2T 
0 (2e - e ) 

-T -2T 
0 (-2e + 2e ) 



/ -T -2T, 

(e - e ) 

-T -2T 
(-e + 2e ) 



A = 



“ -T -2T 

T/2 - 3/4 + e - l/4(e ) 

1/2 - e"^ + l/2(e"^^) 

-T -2T 
e - e 



and, e(kT) = R(kT) (kT) -a^X 2 (kT) -a^X^ (kT) 



The characteristic equation, evaluated at T = 1 second is 
3 2 

Z H- A2^ + A^Z + = 0 , 

with A^ = 0.1998a^ + 0.2325a^ - 1.4191 

A^ = -0.1263a^ - 0.4650a2 + 0.7238 
and A^ = -0.0734a^ + 0.2325a2 - 0.0310 

Observe that it is no longer possible, as was done for the second order 
systems, to include the plant gain, K, in the definition of the system 
parameters. This combination of the gain, K, with the feedback coeffici- 
ents is only possible so long as each coefficient in the numerator of the 
product G(s)H(s) is some product containing a variable feedback gain. In 
the second order case, 

G(s)H(s) = K(a^ + a^s) 

~) 



The characteristic equation can be formed as 

1 + z-transform of [^G(s)H(s)]] 
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Third Order Example 



Thus, the only place that the feedback coefficients occur in the character- 
istic equation, they are each multiplied by the gain, K, However, for 
this third order example, the produce G(s)H(s) has the form, 

G(s)H(s) = K( 1 + a, s + a^s^) 

D (s) 

In the characteristic equation then, K is a multiplier of terms other than 
those involving the feedback coefficients, a^ and a 2 > For the example, we 
are restricted to defining 




Following the usual manipulations, the third order parameter plane, figure 
(15), was achieved. No lines of constant ^ were included on the para- 
meter plane, only to avoid unnecessary cluttering. Observe that the separa- 
tion of complex root regions and real root regions no longer exists; which 
should be expected, as the root locus always has one real root. Also note 
that the region of stability is that region of the parameter plane which 
is enclosed by the contours 1.0, 0^ = -1.0, 0^ = 1.0. From figure 

(15), minimum settling time using feedback compensation is predicted when 

a^ = 1.61 
= 0.54 

The circle of minimum radius enclosing the roots is 0.3. 

Table 5 lists the calculated magnitudes of the state variables at consecu- 
tive sampling intervals in response to a unit step input using these values 
for the feedback coefficients. From Equation (27) then, minimum settling 
time is calculated as 

t = 3T - 4T/ln 0.3 
s 

- 6.3T 
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At time 7T, table 5 shows that the output has settled to within 57o of the 
final value. Figure 16 shows the state variable time response for this 
minimum settling time compensation. These results verify the effective- 
ness of equation (27) in predicting the servo settling time to a good 
engineering approximation, even with the dominance of the roots at maxi- 
mum modulus, as discussed in section 3,4, somewhat uncertain. 

TABLE 5 



nT 


Xj(t) 


X2<t) 


X3(c) 


0 


0 


0 


0 


1 


0.08405 


0.19979 


0.23254 


2 


0.33610 


0.26768 


-0.00653 


3 


0.57736 


0.20644 


-0.06887 


4 


0.74605 


0.13340 


-0.05968 


5 


0.85111 


0.08048 


-0.03964 


6 


0.91357 


0.04724 


-0.02411 


7 


0.95002 


0.02743 


-0.01419 


8 


0 97114 


0.01586 


-0.00824 


9 


0.98334 


0.00916 


-0.00477 


10 


0.99039 


0.00529 


-0.00275 


11 


0.99445 


0.00305 


-0.00159 


12 


0.99680 


0.00176 


-0.00092 


13 


0.99815 


0.00102 


-0.00053 


14 


0.99893 


0.00059 


-0.00031 


15 


0.99938 


0.00034 


-0.00018 


16 


0.99965 


0.00020 


-0.00010 


17 


0.99980 


0.00011 


-0.00006 


18 


0.99988 


0.00007 


-0.00003 


19 


0.99993 


0.00004 


-0.00002 


20 


0.99999 


0.00001 


-0.00001 


A1 gebraic 


Design Methods. 







One simple and quite obvious method for designing feedback compensa- 
tion for minimum prototype response should not be overlooked. In determin- 
ing the system characteristic equation as a function of the feedback vari- 
ables, which must be done to apply the Mitrovic Method, an alternate method 
arises for determining minimum prototype compensation. To obtain minimum 
prototype response for a given system, the function of any compensation 
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Figure 16 

State Variable Time Response for Parameter 
Plane Compensated Third Order System 
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must be to drive all the closed loop roots to the origin of the 2 -plane. 
The characteristic equation of an n-th order system, such as shown on 
figure (10), is 



where 



and 



n 

2 


H- A 2 + 






n-1 


• 0 


A. 

1 


= ®0i^0 + ®li^l •* 


ni 1 



( 28 ) 



n 



a. = feedback coefficient of the i-th derivative 
J 



feedback. Therefore, to algebraically design the minimum prototype compen- 
sation, thereby positioning each closed loop root at the origin of the 2 - 
plane, one has only to simultaneously solve the n non-homogeneous equa- 
tions 

A = 0 
n-1 

\-2 = “ 



Ao " 0 

The solution to these simultaneous equations detexrmines the values of the 
n feedback parameters required to produce minimum prototype response. 
Additionally, as was pointed out in section 3.5, since the number of feed- 
back variables now equals the degree of the characteristic equation, we 
can rewrite the A^ in equation (28) as 

A. = B^.Ka^ + B, .Ka, + + B ,Ka + B , 

1 Oi 0 li 1 ni n n ’ 

thereby providing the designer with yet another potential variable. 

To demonstrate the use of this algebraic design technique, consider 
again the third order example of figure (13). Let the open loop gain be 
variable, K, instead of 1.0 as shown. Furthermore, let the coefficient of 
the position feedback path be a^. Then, for minimum prototype response, 
the characteristic equation becomes 
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where 



1.503 = 0 



3 2 

z + A^z + Aj^z + Aq = 0 

= 0.2325K32 + 0.1998Ka^ + 0.0840K3 q - 
A^ = -0.4650K32 - 0.1263K3^ + 0.1709K3^ + 0.5529 = 0 
Aq = 0.2325K32 - 0.0734K3^ + 0.0188K3^ - 0.0498 = 0 
The simultaneous solution of these equations then, is 

K3q = 3.652 
Ka^ = 4.448 
Ka^ = 1.323 

In order to permit the adoption of the conventional error definition, 
then, K = 3.652 

3o = 1.0 

a^ = 1.218 
= 0.363 

Using these values of the system parameters, the response of the system 
of figure (10) to a unit step input is as shown by table 6. The time 
response of the state variables is shown on figure (17). 



nT 


x^(t) 


TABLE 6 
X2(t) 


X3(t) 


0 


0.000 


0.000 


0.000 


1 


0.307 


0.730 


0.849 


2 


0.929 


0.286 


-0.849 


3 


1.000 


0.000 


0.000 


4 


1.000 


0.000 


0.000 
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Figure 17 

State Variable Time Response for Algebraic 
Compensated Third Order System 
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CHAPTER IV 



DIGITAL CCMPUTATION IN MULTI -ORDERED PARAMETER SPACE 

While the feedback compensation of many n-th order sampled-data con- 
trol systems can be reduced to the two parameter problem presented in 
Chapter III, circumstances can, and do, arise where derivative signals 
higher than the second derivative are available. Since each additional 
degree of freedom available to the designer enlarges the region in the z- 
plane in which the closed loop roots may be positioned, common-sense alone 
dictates their usage where available. As pointed out in preceding chapters, 
Mitrovic*s Method is primarily two parameter in nature. Additional design 
techniques are then desirable whereby a parameter space higher than a 2- 
space can be achieved. 

The digital computer provides one method of increasing the order of 
the parameter space. Any design process or technique is simply a method 
for locating the roots of a polynomial, the system characteristic equation, 
as a function of one or more parameters. These methods were developed pri- 
marily to surmount the difficulty of factoring a given polynomial. With 
the digital computer however, this polynomial factoring can be repeatedly 
performed, always with greater accuracy, and sometimes in less total time 
than analysis performed by the graphical and semi-graphical techniques. In 
this chapter, one approach to multi-parameter study using the digital com- 
puter will be presented. The emphasis on the application of this method 
will be directed toward feedback compensation; however, the extension to 
cascade compensation or other parameter variations does exist. 

Recent emphasis on the state space method of analysis, [2], as applied 
to sampled-data systems affords an excellent, if not necessary, point of 
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departure for digital computer analysis. The general form of the state 
space vector-matrix equation describing a sampled-data feedback control 
system was given by equation (14), 

[(k + 1)T] = $ X(kT) + A R(kT) 

It can easily be demonstrated, see appendix 1, that the characteristic 
equation of the state transition matrix, is indeed the characteristic 

equation of the closed-loop transfer function. Since numerous numerical 
methods presently exist for determining the eigenvalues of a given matrix, 
the determination of the closed loop pole locations by evaluating the 
eigenvalues of the ^ matrix is readily performed on the digital computer. 
To determine these root locations as a function of one or more parameters, 
one has only to repeatedly increment the parameters and determine the eigen- 
values of the resulting ^ matrix. In this fashion, an n-coordinate 
space can readily be created by computer programming techniques. One 
method is to use the FORTRAN language, and program DO LOOP'S inside other 
DO loop's until sufficient coordinates are defined to permit the designer 
to sweep through a desired region in the n-space, determining the root 
locations corresponding to each point in the n dimensional parameter space. 

In order to demonstrate this method of analysis more fully, the 
following specific project was undertaken. Given the perfectly general 
second order sampled-data system shown in figure (18) , the values of the 
parameters Ka^ and Ka^ required for minimum prototype, minimum settling time 
response to a step input are to be determined as functions of the sampling 
interval, T. The resulting computer program should be such that the de- 
sired design data can be readily generated with only a minimum amount of 
additional input data supplied to the data, for any second order system 
with real poles. 
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The system of figure (18) can be represented by the open-loop vector 
matrix equation: 

X[(k 1)T] = FX(kT) -f A e(kT). 

The error signal at the k-th sampling interval is represented as; 

e(kT) = R(kT) - a^X^(kT) - a^X^CkT) 

In the chosen example, the F matrix is perfectly general, irregardless of 
the location of the olant poles and zeros, and 



-p,T "p„T 
pel -p, e 2 
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e 



- e 



■P 2 T 



’" 2 ^'^ - 1)/0q 



•P^e'^l^ + P2 s'^2^ 



where 



Modifications to the F matrix are needed, however, if a pure second order 
system is used. The ^ matrix on the other hand, is dependent upon the 
number of pure integrations in the forward path of the plant transfer 
function. For example, for ~ ^2 ~ 



A =K 



+ b^T bpT^/2 



i 



(T) + b, + b^T 



1 



0 



When only one pure integration exists, that is p^^O, p^ = 0, 

(1/pp^ [(bj + b^T)p^ - b^ + O^e’Pl^] 



Arr K 



where , 



S’ (I) 



+ b^/p^ - 



^2 



(b2Pi - bpp^ + b^ 
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Block Diagram 
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General Second Order 
System With Zeros 
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Finally, when n^ ^ 0, and P 2 ^ 



K 



where 






b2S^(T) + G^Q^e'Pl^ -G20^e‘P2'^ 



A 2 

Q 3 = l/(Pj - PjP2^ 



Gj S pj - Pjbj + b(, 



°2 ' »2 - Pa\ ’’ 



0 



Lse. 



0 ^ = 1/(P2 - P 1 P 2 ) 

^5 ^ ^^2’^2 ■ ’^1^^2 '*' impuli 

The program which was developed, see appendix 2, to satisfy the 

design goal requires four inputs in the form of data. In the order in 
which they are required, they are: 

a. the pole locations of the plant in the s-plane. The nega- 
tive value of this s-plane location must be entered on the data card, and 
the poles should be listed in decreasing magnitude order. 

b. the plant s-plane numerator coefficients listed in ascend- 



ing order. 

c. the minimum and maximum values of T, the sampling interval, 
and the incremental value of the sampling interval to be used in the in- 
vestigation. T minimum must be greater than zero. 

d. miscellaneous title information and graph scales required 
for graph output in accordance with an existing graph subroutine (DRAW) . 

The program outputs are: 

a, a tabulated listing of the minimum settling time values of 
Ka^ and Ka^ determined for each value of T used in the investigation. 

b. a graph plot of the minimum settling time value of Ka^ vs. 



T. 



c. a graph plot of the minimum settling time value of Ka^ vs. 



T. 
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Let us now consider the geometric significance of this program. It 
has been demonstrated that for a given value of the sampling interval T, 
the origin of the ^-plane maps into a single point on the parameter plane 
for a second order system. However, if a three dimensional right-handed 
coordinate set is created, whose axes are Ka^^ , ^^ 2 * minimum 

settling time point corresponding to = 0 will map into a curve in the 
3-space. The plots of Ka^ and Ka^ required for minimum settling time re- 
sponse vs. T can be thought of as the projection of this curve in 3-space 
onto the Ka^^-T and Ka^-T planes. Similiarly, any n-th order system, with 
m degrees of freedom, where m ^ n, can be investigated in a computer 
oriented m-space and displayed to the designer by a set of m projections 
on the a^ vs. T, (i = 1, 2,...,m), planes. If m = n, a set of projections 
on the Ka^-T planes, (i = 1, 2,...,n), are available. 

The following examples will serve to illustrate the class of feed- 
back compensation problems, the solutions of which can be greatly simpli- 
fied by using this particular programming technique. 

EXAMPLE #1 

A sampled-data control system whose plant transfer function is 

G(s) = K(s + l)/s(s 4- 2) , 

is to be used with a dual capacity for sampling rates of lOcps and Icps. 
Compensation is to be designed to provide minimum prototype response to 
a step input for each of the possible sampling frequencies. 

Figures (19) and (20) are the minimum settling time parameter planes 
for this plant transfer function. From these parameter planes, the re- 
quired parameter values for each sampling frequency can be obtained as: 

T Ka. ^ 

0.1 sec 10.0^ 0.949 

1.0 sec 1.15 0.538 
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Finally, for unity position feedback, the physical plant parameters 



must be: 



T ^2 

0.1 sec 10.0 0.095 

1.0 sec 1.15 0.467 

Figure (21) represents the block diagram of the compensated system. 

When the three switches in the physical plant are ’’ganged”, then the control- 
led switch of the sampling frequency will introduce the correct, associated 
feedback compensation for minimum settling time. 
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EXAMPLE #2 



A system whose plant transfer function is 

G(s) = K(s + Z)/s(s + 2) 

is to be controlled by a sampled error signal. The plant zero location 
is a function of a cascaded compensator, and can be placed an)rwhere on 
the real axis of the s-plane. A sampling interval of 0.5 seconds is 
intended. The designer is to set the gain, K, the zero location, Z, and 
design feedback compensation for minimum settling time in response to a 
step input, and minimum steady state error for a ramp input. 

Using the program to determine the minimum settling time values of 
Ka^ and Ka 2 as a function of T, a family of curves, with parameter Z, 
can be constructed on the Ka^ vs. T and Ka^ vs. T planes. These families 
of curves are shown by figures (22) and (23). Of these two families of 
curves, the minimum settling time Ka^ vs. T plane is by far the most 
interesting. The minimum settling time Ka^^ vs. T loci was shown by fig- 
ure (19) to be hyperbolic. Furthermore, section 3.3 outlined the presence 
of discontinuities on parameter planes for systems with more than two 
parameters. These discontinuities again appear on the minimum settling time 
parameter planes, figures (22) and (23), and are best demonstrated by fig- 
ure (23). Inspection of figure (22), the minimum settling time Ka^ vs. T 
plane, also reveals these discontinuities on the loci of minimum settling 
time for the cases of Z less than zero (right-half s-plane zeros). 

It has been demonstrated, [2] and [7], that for a type one sampled- 
data control system, the steady state error in response to a ramp input is 
inversely proportional to the plant root locus gain; similiar to the *'velo- 
city lag** error in continuous systems. Consequently, to minimize the ramp 
error and at the same time design for minimum settling time response to a 
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step input, the designer must maximize the parameter K (or Ka^ when using 
unity position feedback) on the specified line of constant sampling in- 
terval T = 0,5. Inspection of figure (22) shows the importance played in 
this selection of the plant gain by the discontinuity phenomena. From 
figure (22), the maximum Ka^^ attainable for T = 0.5 is Ka^ = 12.7, when 
Z =-1.5. For Z less than -1.5, the discontinuity asymptote lies to the 
left of T = 0.5, and a negative value of Ka^ (positive feedback) is re- 
quired for minimum step input settling time. For Z greater than -1.5, 
lower values of Ka^ are required for minimum settling time for a step in- 
put and therefore a larger ramp error is experienced. 

Therefore, the design criteria of this example is achieved for: 

K = 12.7 
a^ = 1.0 
a2 = 0.656 
Z = -1.5 
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CHAPTER V 



CONCLUSIONS 

The application of Mitrovic's Method to feedback compensation of 
sampled-data control systems has been shown to be quite useful in the 
minimization of system settling time* The parameter plane for the second 
order system presents an easily obtainable picture of the system stability 
and settling time as a function of the system parameters- Using the 
methods described herein, the second order parameter plane development 
results in a two parameter analysis method requiring about the same total 
effort to construct as conventional one parameter methods* 

In the case of higher order systems, the effectiveness of the Mitrovic 
Method as a design tool is somewhat reduced* Foremost, the plotting of the 
various loci requires at least a desk calculator, preferably a digital 
computer. If the number of system parameters can be reduced to only two, 
then the parameter plane is still a good tool for minimizing settling time* 
Unfortunately, the determination of the minimum settling time compensation 
from the higher order parameter plane requires a more tedious study of the 
plane than for the second order case. The intersections of all the loci 
of constant and constant ^ must be examined in detail in search of 

the circle of minimum radius in the z-plane encircling all the closed loop 
roots . 

For higher order systems with more than two degrees of freedom, 
analysis with the digital computer is most effective, indeed necessary* 

The vector-matrix state space equations representing the given system 
were shovm to lead directly to the determination of the closed loop root 
positions as a function of up to n parameters, n being the order of the 
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system. The strength of this method is the ability to reduce the design 
criteria to a function of almost any physical parameter or condition. 
However, since a large fraction of the feedback compensation problems are 
not more than two parameter in nature, the Mitrovic Method can be an 
exceedingly important tool to the design engineer. 
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APPENDIX I 



STABILITY ANALYSIS OF VECTOR -MATRIX DIFFERENCE EQUATIONS 
The general form of the state space vector-matrix equation describ- 
ing a linear system with a sampled input is, 

2^[(k + 1)T] = FMkT) -h A u(kT). (I-l) 

In this form, the states of the system, X, are defined as the control 
system output and its first n-1 derivatives, where n is the order of the 
linear system. For a feedback control system, the control vector, u(t), 
is the difference between the command input, R(t), and the sum of the 
feedback quantities, ie; 

u(t) = R(t) - X(t), (1-2) 

where A*" is the transpose of A, the coefficient vector controlling the 
magnitude of feedback for each state variable. In discrete time form, 
equation (1-2) is written as 

u(kT) = R(kT) - X(kT) . (1-3) 

By substituting equation (1-3) into equation (I-l), the familiar vector- 
matrix difference equation for a sampled data control system is obtained 
as 

X[(k + 1)T] = ^'^((kD-AA*^ X(kT) + AR(kT) , (1-4) 

or by imbeding the feedback quantities directly into the F matrix, 

X[(k + 1)T] = J X(kT) + A R(kT) . (1-5) 

With the feedback thus imbeded in the transition matrix, ^ , the 
Stability analysis of a given system can be confined to an analysis of the 
partial state equation, 

X[(k + 1)T] X(kT) , (1-6) 

which is an initial value problem. If the state trajectory returns to a 
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givon erjuilibrium, or singular point, in the state space for any initial 
stc'^te vector, it will also move to a predictable translation of this singu- 
lar point and achieve equilibrium for the input R(t) which can be described 
by the initial state vector. The dynamic stability of the system is thus 
ensured . 



Hahn, [1], studied the stability of equation (1-6) using Liapunov’s 
Direct Method, and concluded that the stability of equation (1-^6) was 
ensured when the logarithm of the eigenvalues of the J matrix were all 
less than zero in magnitude, that is 

In Xt < 0 • (1-7) 

In (Re X|_+ j ImXl.) < 0 

In ( I X(( ) <0 

In I + j ( 0 ± 2n7r)<0 

which is satisfied if and only if 

< 1 

Immediately, a direct analogy of the unit circle in the - plane can be 

drawn to the unit circle of the z-plane. This analogy can also be demon- 
strated by considering the z-t ransformat ion of equation (1-6) as: 

zX(z) = ^ X(z) -I- zX(0) 

(zl - 3 )X(z) = zX(0) 

or, ^(2^) z(zl - ^ ^ 

Now, the inverse of the matrix (zl -3^) is the transpose of the cofactor 
matrix, divided by the determinant value of (zl - 5 ); or, if the cofactor 
transpose matrix is denoted as [C] 



This is equivalent to 

or 

or 



(zl - ^ ) 



[£2 

1 2l - £ I 
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Then 



X(^) = 



z[C] 



X(0) 



|zl - J I 

Since z is just a scalar variable over a complex region, the determinant 
( zl * ^ I is nothing more than the characteristic equation of the 
matrix, whose characteristic values or eigenvalues are the closed loop 
root locations in the z-plane, or ^ -plane. Similiarly, in equation 
(I-l), the eigenvalues of the F matrix are the open loop root locations 
in the z-plane ( \ -plane). The vector-matrix formulation is ideally 
suited for digital computation of the state variable time response in 
discrete form. Even more important, it is also suited for determination 
of the system dynamic stability by digitally locating the eigenvalues of 
the J matrix as a function of any variable in the ^ matrix. 
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APPENDIX II 



FORTRAN PROGRAM FOR MINIMUM SETTLING TIME PARAMETER PLANE 



PROGRAM design 

C PROGR^^'MER L.D.NACF 

C TH-IS PROGRAM COMPUTES AND GRAPHS THE VALUES OF THE 

C PARAMETERS KAl AND KA2 REQUIRED FOR MINIMUM PROTOTYPE 

C RESPONSE OF ANY SECOND ORDER SAMPLED-DATA FEEDBACK 

C CONTROL SYSTEM AS A FUNCTION OF THE SAMPLING INTERVAL T. 

C INPUT DATA IS 

C CARD 1 VALUES OF P» THE NEGATIVE OF THE PLANT S-PLANE 

C POLE LOCATION IN DECREASING MAGNITUDE ORDER* FORMAT 2FI0.5 

C CARO 2 COEFFICIENTS OF THE S-PLANE NUMERATOR OF PLANT 

C IN ASCENDING ORDER* FORMAT 3F10.5 

C CARD 3 T MIN* T MAX* DELTA T* IN FORMAT 3F10.5. T MIN IS 

C THE MINIMUM VALUE OF T SCANNED. T MAX IS THE MAXIMUM 

C VALUE. AND DELTA T IS THE INCREMENTAL VALUE OF T USED IN 

C cA^LCULATIONS. MAXIMUM NUMBER OF POINTS COMPUTED IS 500. 

C T MIN MUST BE GREATER THAN ZERO. 

C CARD A FIRST LINE OF THE GRAPH TITLE* COLUMNS 1-48. 

C CARD 5 SECOND LINE OF THE GRAPH TITLE* COLUMNS 1-48 

C CARD 6 EXSCALE.YSCALE FOR GRAPH OUTPUT IN FORMAT 2F10.5 

DIMENSION P{ 2) *0( 3 ) * ANS 1 ( 500 ) * ANS2 ( 500 ) * T T ( 500 ) 

DIMENSION ITITLEI12) ^ 

READ lOO.P o 

100 FORMAT (2F10.5) 

75 FORMAT (/////) 

76 FORMAT (//) 

PRINT 75 
PRINT 101 

101 FORMAT (47HVALUES OF P, NEGATIVE OF S-PLANE POLE LOCATIONS) 
PRINT 76 

PRINT 100, ( P( I ) , 1 = 1 ,2 ) 

READ 102*8 

102 FORMAT (3F10.^) 1 

PRINT 75 
PRINT 103 

103 FORMAT (47HPLANT NUMERATOR COEFFICIENTS IN ASCENDING ORDER) 
PRINT 76 

PRINT 102*(5(I) *1=1*3) 

READ 102* TMIN.TMAX.DELTAT 
PRINT 104 

104 FORMAT (/////* 3X * 5HT MIN*6X*5HT MAX ♦ 6X * 7HDELTA T*//) 

PRINT 102* TMIN*TMAX*DELTAT 

200 FORMAT (6A8) 

READ 200* ( ITITLEI I ) * 1 = 1*6) j 

READ 200* ( ITITLEI I )* 1 = 7*12 ) v 

PRINT 201 

201 FORMAT (/////* 12HGRAPH TITLES*//) 
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PRINT 200 ♦ ( I TI TLE ( I ) » I = 1 »6 ) 

PRINT 20n,( ITITLEI 1 ) »I=7»12) 

READ 202» EXSCALEjYSCALE 

202 FORMAT (2F10.5) 

PRINT-203 

203 FORMAT (/////» lOX ♦ 7HX-SCALE » lOX ♦ 7HY-5CALE »// ) 
PRINT 20A, EXSCALE »YSCALE 

20A FORMAT ( 8 X , FI 0 . 5 ♦ 9X ♦ F 1 0 . 5 ) 

LAST = 0 ■ 

P1=P(1) 

P2=P(2) 

B0=3 ( 1 ) 

B1=B (2 ) 

B2=B(3) 

K = 0 . • 



> I 

IB; 



I' 

b.u 
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fe 



105 



118 



I 



106 

107 

108 

109 

110 



111 



T=TMIN 
K = K + 1 

IF (P1+P2) 118»112»118 

U=1 . /EXPF (P1*T) 

V=1./EXPF(P2*T) 

PHI 1=P2^^U/ ( P2-P1 ) +P1*V/ (P1-P2 ) 

PH12=U/(P2-P1 )+V/(Pl-P2) 

PH2‘l = ( U«V-1 , )*( (P1*P2 ) / (P1 + P2 ) ) 

PH2 2 = P1*-U/ ( P1-P2) +P2*V / (P2-P1 ) 

N = C 

IF (PI) 107,106»107 
N = N+1 

IF (P2) 109, 108, 109 

N = N+1 

IF (N-1) 110,111,112 

D1 = B0/ ( P1*P2 ) + ( (B2*P1*P1-81^^P1 + B0) / ( P1*P1-P1*P2 ) ) *U 
l+( ( B2*P2*P2-B1*P2+B0 ) / (P2^P2-P1*P2 ) )*V 
D2= ( ( B1-P1-P2 )■» (-P1 ) +B0-P1*P2 )*U/ ( P2-P1 ) 
l+( (B1-P1-P2 )*(-P2)+30-Pl*P2)*V/(Pl“P2 ) 

GO TO 113 

D1 = ( (B1+B0*T)*P1-B0) / (P1*P1 )+( ( B2*P1*P1“B1*P1+B0)*U) / (P1*P1 ) 
D2=B0/Pl-( ( B2*P1-B1 ) *P1+B0 ) *U/P1 
GO TO 113 



112 Dl=B 2 + Bl^^T + B 0 *T*T/ 2 . 
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llA 



D 2 =B 1 +B 0 *T I 

PH 11 = 1 . 

PH 12 =T 

PH 21 = 0 . 

PH 22 = 1 . 

DEN 0 M=D 1 * ( PH 21 *D 1 -PH 11 *D 2 )-D 2 *( PH 12 *D 2 -PH 22 *D 1 ) 

DUMA 1 = ( PHI 1 + PH 2 2 ) * (PH 21 *D 1 -PH 1 1 *D 2 ) 

1 +D 2 ^f (PH 11 *PH 22 -PH 12 *PH 21 ) 

DUMA 2 = -D 1 *( PHI 1 *PH 22 -PH 12 *PH 21 )'-(PH 11 + PH 22 )*(PH 12 *D 2 -PH 2 2 *D 1 ) 
ANSI (K) =DUMA 1 /DEN 0 M 
ANS 2 (K)=DUMA 2 /DEN 0 M 
TT(K)=T 

IF (TMAX-T) 115 , 115,114 

T=T+DELTAT 

GO TO 105 
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115 PRINT 116 

116 FORMAT (/////»10X»1HT»19X»3HI^A1,17X»3HKA2»// ) 

PRINT 117» (TT( J) »ANS1 (J ) .ANS2I J) »J = 1»K) 

117 FORMAT! 5X,F 10.5 >10X, FI 0.5 »10X, Flo. 5) 

CALL DRAW ( K * T T . ANSI » 0 » 0 ♦ AHKAl » I T I TLE » EXSCALE » YSCALE . 
1G,0,0,0»8,8»1»LAST) 

CALL DRAW, (K»TT»ANS2»0»0»4HKA2 . I T I TLE » EXSCALE » YSCAL E ♦ 
ID *'0,0, 0,8, 8,1 .LAST) 

END 

END 
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